#!/usr/bin/python3
import sys
import argparse
import copy
import hisdheader
import struct
import netCDF4 as nc 
from hisd_pixlin2lonlat import HisdPixLin2LonLat
from getdata_by_pixlin import GetData
from myerror import *
import time
import numpy

MAXFILE=10
INVALID=-1

WIDTH=5501	#/* default pixel number */
HEIGHT=2001	#/* default line number */
LTLON=90.0	#/* default left top longitude */
LTLAT=10.0	#/* default left top latitude */
DLON=0.01	#/* default Spatial resolution (longitude) */
DLAT=0.01	#/* default Spatial resolution (latitude) */


class Hisd2netcdf:
    def __init__(self,args) -> None:
        self.param={
            "ltlat":args.lat[0], #left top longitude
            "ltlon":args.lon[0], #left top longitude
            "width":args.width[0], #pixel number
            "height":args.height[0],#line number
            "dlon":args.dlon[0],    #Spatial resolution (longitude)
            "dlat":args.dlat[0]     #Spatial resolution (latitude)
        }
        self.arg={
            "outfile":args.o[0], #output file 
            "infile":args.i,
            "filenum":len(args.i) #input file
        }
        self.data={}
        self.data["phys"]=[]
        if args.width[0]<10:
            self.param["width"]=10
        if args.height[0]<10:
            self.param["height"]=10
        if args.lat[0]<-90.0 or 90.0<args.lat[0]:
            self.param["ltlat"]=LTLAT
        if args.lon[0]<-180.0 or 180.0<args.lon[0]:
            self.param["ltlon"]=LTLON
        if args.dlat[0]<0.0 or 10.0<args.dlat[0]:
            self.param["dlat"]=DLAT
        if args.dlon[0]<0.0 or 10.0<args.dlon[0]:
            self.param["dlon"]=DLON

        
        
    def getData(self):
        n=-1
        minLine=99999.0
        maxLine=-99999.0
        fp=[]
        startLine=[]
        endLine=[]
        hisd=[]
        n_size=self.param["width"]*self.param["height"]
        getdata=[]
        pixlin2lonlot=[]
        #初始化一个空列表
        for i in range(n_size):
            pixlin2lonlot.append(None)
        for i in range(self.arg["filenum"]):
            #1-1 open file
            try:
                fp_temp=open(self.arg["infile"][i],"rb")
            except Exception:
                sys.stderr.write("cannot open file.")
                continue
            fp.append(fp_temp)
            hisd.append(hisdheader.HisdHeader(fp[i]))
            #1-2 read hisd header
            if NORMAL_END!=hisd[i].read():
                error="error : read header ["+self.arg["infile"][i]+"]\n"
                sys.stderr.write(error)
                continue
            #1-3 starLine and endLine
            startLine.append(struct.unpack("H",hisd[i].segm_info["strLineNo"])[0])
            endLine.append(startLine[i]+int.from_bytes(hisd[i].data_info["nLin"],sys.byteorder)-1)
            #1-4 check header consistency
            if n==-1:
                n=i
            if hisd[n].calib_info["bandNo"]!=hisd[i].calib_info["bandNo"] or hisd[n].calib_info["gain_cnt2rad"]!=hisd[i].calib_info["gain_cnt2rad"] or hisd[n].proj_info["loff"]!=hisd[i].proj_info["loff"] or hisd[n].proj_info["coff"]!=hisd[i].proj_info["coff"]:
                sys.stderr.write("header consistency error\n")
                error=self.arg["infile"][n]+":"+self.arg["infile"][i]+"\n"
                sys.stderr.write(error)
                return ERROR_INFO
            n=i
        # 1-5 check file open
        if n==-1:
            sys.stderr.write("error : can not open all files\n")
            for f in fp:
                f.close()
            return ERROR_FILE_OPEN
        #1-6 satellite name & band number 
        self.param["band"]=struct.unpack("H",hisd[n].calib_info["bandNo"])[0]
        self.data["satName"]=hisd[n].basic_info["satName"].decode()
        for l in range(self.arg["filenum"]):
            getdata.append(None)
        #2 get data
        for j in range(self.param["height"]):
            for i in range(self.param["width"]):
                #2-1 init
                count=struct.unpack("H",hisd[n].calib_info["outCount"])[0]
                k=j*self.param["width"]+i
                #2-2 convert lon & lat to pix & lin 
                pixlin2lonlot[k]=HisdPixLin2LonLat(hisd[n])
                pixlin2lonlot[k].lonlat_to_pixlin(self.data["lon"][i],self.data["lat"][j])
                #2-3 min & max line
                if pixlin2lonlot[k].pix==-9999 or pixlin2lonlot[k].lin==-9999:
                    self.data["phys"][k]=INVALID
                    continue
                if minLine>pixlin2lonlot[k].lin:
                    minLine=pixlin2lonlot[k].lin
                if maxLine<pixlin2lonlot[k].lin:
                    maxLine=pixlin2lonlot[k].lin
                #2-4 get count value
                hisd[n].correcttable["flag"]=0
                    
                for l in range(self.arg["filenum"]):
                    if startLine[l]-0.5<=pixlin2lonlot[k].lin and pixlin2lonlot[k].lin<endLine[l]+0.5:
                        getdata[l]=GetData(hisd[l])
                        getdata[l].getdata_by_pixlin(pixlin2lonlot[k].pix,pixlin2lonlot[k].lin)
                        count=getdata[l].sout
                        break
                #2-5 check count value
                if count==int.from_bytes(hisd[n].calib_info["outCount"],sys.byteorder) or count==int.from_bytes(hisd[n].calib_info["errorCount"],sys.byteorder):
                    self.data["phys"][k]=INVALID
                else:
                #2-6 convert count value to radiance
                    #print(count)
                    radiance=float(count)*struct.unpack("d",hisd[n].calib_info["gain_cnt2rad"])[0]+struct.unpack("d",hisd[n].calib_info["cnst_cnt2rad"])[0]
                    #2-6 convert radiance to physical value
                    if (struct.unpack("H",hisd[n].calib_info["bandNo"])[0]>=7 and hisd[n].basic_info["satName"].decode().find("Himawari")!=-1) or (struct.unpack("H",hisd[n].calib_info["bandNo"])[0]>=2 and hisd[n].basic_info["satName"].decode().find("MTSAT-2")!=-1):
                        phys=hisd[n].radiance_to_tbb(radiance)
                        self.data["phys"][k]=phys
                    else:
                        # visible or near infrared band 
                        self.data["phys"][k]=struct.unpack("d",hisd[n].calib_info["rad2albedo"])[0]*radiance
        #3 convert maxLine & minLine to scanTime 
        for l in range(self.arg["filenum"]):
            #3-1 startTime
            if startLine[l]<=minLine and minLine<=endLine[l]:
                for i in range(1,struct.unpack("H",hisd[l].obstime["obsNum"])[0]):
                    if minLine<struct.unpack("H",hisd[l].obstime["lineNo"][i])[0]:
                        self.data["startTime"]=struct.unpack("d",hisd[l].obstime["obsMJD"][i-1])[0]
                        break
                    elif minLine==struct.unpack("H",hisd[l].obstime["lineNo"][i])[0]:
                        self.data["startTime"]=struct.unpack("d",hisd[l].obstime["obsMJD"][i])[0]
                        break
            #3-2 endTime
            self.data["endTime"]=0
            if startLine[l]<=maxLine and maxLine<=endLine[l]:
                for i in range(1,struct.unpack("H",hisd[l].obstime["obsNum"])[0]):
                    if maxLine<struct.unpack("H",hisd[l].obstime["lineNo"][i])[0]:
                        self.data["endTime"]=struct.unpack("d",hisd[l].obstime["obsMJD"][i-1])[0]
                    elif maxLine==struct.unpack("H",hisd[l].obstime["lineNo"][i])[0]:
                        self.data["endTime"]=struct.unpack("d",hisd[l].obstime["obsMJD"][i])[0]
                        
        #4 check data
        print("Satellite Name:",self.data["satName"],sep="")
        print("Band Number    :",self.param["band"],sep="")
        print("physical value :")
        j=0
        while j<self.param["width"]:
            print("%6.1f "%self.data["lon"][int(j)],end=" ")
            j=j+(self.param["width"]-1)/20
        print("\n",end="")
        i=0
        while i<self.param["width"]:
            k=i*self.param["width"]+j
            print("%6.1f "%self.data["lat"][int(i)],end=" ")
            j=0
            while j<self.param["width"]:
                k=i*self.param["width"]+j
                print("%6.2f "%self.data["phys"][int(k)],end=" ")
                j=j+(self.param["width"]-1)/20
            print("\n",end="")
            i=i+(self.param["width"]-1)/20
        # 6 fclose
        for f in fp:
            f.close()
        return NORMAL_END
    def __defNetCDF(self,nf) -> int:
        ch=self.param["band"]
        satName=self.data["satName"]
        outFile=self.arg["outfile"]
        physDimID=[]
        latName="latitude"
        lonName="longitude"
        startTimeName="start_time"
        endTimeName="end_time"
        timeUnit="days since 1858-11-17 0:0:0"
        invalid=float(INVALID)
        self.ncparam={}
        # 1 phys name 
        if ch<6:
            physName="albedo"
            physUnit="1"
            physStdName="reflectivity"
            tittle=satName+" band-"+str(ch)+" ALBEDO"
        elif ch <= 16:
            physName="tbb"
            physUnit="K"
            physStdName="brightness_temperature"
            tittle=satName+" band-"+str(ch)+" TBB"
        else:
            sys.stderr.write("ch number error  (ch number : 1 to 16)\n")
            return ERROR_INFO
        # 2 dimensions 
        latDimIns=nf.createDimension(latName,self.param["height"])
        lonDimIns=nf.createDimension(lonName,self.param["width"])
        startTimeDimIns=nf.createDimension(startTimeName,1)
        endTimeDimIns=nf.createDimension(endTimeName,1)
        # 3 variables
        self.ncparam["latVarIns"]=nf.createVariable(latName,"f4",(latName,)) 
        self.ncparam["lonVarIns"]=nf.createVariable(lonName,"f4",(lonName,)) 
        self.ncparam["physDimIns"]=[latDimIns,lonDimIns]
        self.ncparam["physVarIns"]=nf.createVariable(physName,"f4",(self.ncparam["physDimIns"][0],self.ncparam["physDimIns"][1]),fill_value=invalid)
        self.ncparam["startTimeVarIns"]=nf.createVariable(startTimeName,"f8",())
        self.ncparam["endTimeVarIns"]=nf.createVariable(endTimeName,"f8",())
        # 4 attribute
        self.ncparam["latVarIns"].units="degrees_north"
        self.ncparam["latVarIns"].long_name="latitude"
        self.ncparam["lonVarIns"].units="degrees_east"
        self.ncparam["lonVarIns"].long_name="longitude"
        self.ncparam["physVarIns"].units=physUnit
        self.ncparam["physVarIns"].long_name=physStdName
        self.ncparam["startTimeVarIns"].units=timeUnit
        self.ncparam["startTimeVarIns"].standard_name="time"
        timeName="observation start time"
        self.ncparam["startTimeVarIns"].long_name=timeName
        self.ncparam["endTimeVarIns"].units=timeUnit
        self.ncparam["endTimeVarIns"].standard_name="time"
        self.ncparam["endTimeVarIns"].long_name="time"
        timeName="observation end time"
        self.ncparam["endTimeVarIns"].long_name=timeName
        # 5 global attributes
        nf.tittle=tittle
        nf.institution="MSC/JMA"
        cc=satName+" satellite observation"
        nf.source=cc
        # 6 hh:nn:ss mm/dd/yyyy: created.
        num=outFile.rfind("/")
        if num==-1:
            fileName=outFile
        else:
            fileName=outFile[num:]
        iDate=time.gmtime()
        strLength=38+len(fileName)+9+len(nc.getlibversion())+1
        str1=time.strftime("at %H:%M:%S %m/%d/%Y: file created. ",iDate)
        history=str1+fileName+nc.getlibversion()
        nf.history=history
        nf.Conventions="CF-1.4"
        # 7 Quit 
        return NORMAL_END
    def __putNetCDF(self,nf):
        self.ncparam["latVarIns"][:]=self.data["lat"]
        self.ncparam["lonVarIns"][:]=self.data["lon"]
        phys=numpy.zeros((self.param["height"],self.param["width"]),dtype="float32")
        for i in range(self.param["height"]):
            for j in range(self.param["width"]):
                phys[i][j]=self.data["phys"][i*self.param["height"]+j]
        self.ncparam["physVarIns"][:]=phys
        self.ncparam["startTimeVarIns"][:]=self.data["startTime"]
        self.ncparam["endTimeVarIns"][:]=self.data["endTime"]
        return NORMAL_END
    def makeNetCDF(self):
        # 1 open netcdf file 
        nf=nc.Dataset(self.arg["outfile"],"w")
        # 2 define netcdf file 
        self.__defNetCDF(nf)
        # 3 put data
        self.__putNetCDF(nf)
        # 4 close netcdf file 
        nf.close()
        # 5 end 
        return NORMAL_END
        
        

# 1 get arguments 
parser=argparse.ArgumentParser()
parser.add_argument("-o",nargs="*",default=None)
parser.add_argument("-i",nargs="*",default=None)
parser.add_argument("-width",nargs="*",default=WIDTH,type=int)
parser.add_argument("-height",nargs="*",default=HEIGHT,type=int)
parser.add_argument("-lat",nargs="*",default=None,type=float)
parser.add_argument("-lon",nargs="*",default=None,type=float)
parser.add_argument("-dlat",nargs="*",default=DLAT,type=float)
parser.add_argument("-dlon",nargs="*",default=DLON,type=float)
args = parser.parse_args()

# 2 check error   
if hasattr(args,"o")==False or hasattr(args,"i")==False:
    cPtrnum=sys.argv[0].rfind('/')
    cPtr=sys.argv[0][cPtrnum:]
    error="Usage:"+cPtr+'''
            [OPTION]\n
			  -i <InFile> [-i <InFile2> ...]\n
			  -o <OutFile>\n
			  -width  <Pixel Number>\n
			  -height <Line Number>\n
			  -lat    <Left top latitude>\n
			  -lon    <Left top longitude>\n"
			  -dlat   <Spatial resolution (longitude)>\n
			  -dlon   <Spatial resolution (latitude)>\n
            '''
    sys.stderr.write(error)
h2n=Hisd2netcdf(args)
# 3 check parameter
print("Left top (lat,lon) : (%6.2f,%6.2f)\n"%(h2n.param["ltlat"],h2n.param["ltlon"]))
print("width,height       : (%6d,%6d)\n"%(h2n.param["width"],h2n.param["height"]))
print("Spatial resolution : (%6.2f,%6.2f)\n"%(h2n.param["dlat"],h2n.param["dlon"]))
n_size=h2n.param["width"]*h2n.param["height"]

# 4 init
h2n.data["phys"]=[]
h2n.data["lat"]=[]
h2n.data["lon"]=[]
for i in range(h2n.param["height"]):
    h2n.data["lat"].append(0.0)
for i in range(h2n.param["width"]):
    h2n.data["lon"].append(0.0)
for i in range(n_size):
    h2n.data["phys"].append(INVALID)
# 5 set longitude and latitude
for i in range(h2n.param["height"]):
    h2n.data["lat"][i]=h2n.param["ltlat"]-h2n.param["dlat"]*i
for i in range(h2n.param["width"]):
    h2n.data["lon"][i]=h2n.param["ltlon"]+h2n.param["dlon"]*i
# 6 get data
h2n.getData()
# 7 make NetCDF
h2n.makeNetCDF()
# 8 End
print("Normal End")